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ABSTRACT 

Analysis of high-precision timing observations of an array of '^20 millisecond pulsars (a so- 
called "timing array") may ultimately result in the detection of a stochastic gravitational-wave 
background. The feasibiUty of such a detection and the required duration of this type of ex- 
periment are determined by the achievable rms of the timing residuals and the timing stability 
of the pulsars involved. We present results of the first long-term, high-precision timing cam- 
paign on a large sample of millisecond pulsars used in gravitational-wave detection projects. 
We show that the timing residuals of most pulsars in our sample do not contain significant low- 
frequency noise that could limit the use of these pulsars for decade-long gravitational-wave 
detection efforts. For our most precisely timed pulsars, intrinsic instabilities of the pulsars or 
the observing system are shown to contribute to timing irregularities on a five-year timescale 
below the 100 ns level. Based on those results, realistic sensitivity curves for planned and 
ongoing timing array efforts are determined. We conclude that prospects for detection of 
a gravitational-wave background through pulsar timing array efforts within five years to a 
decade are good. 
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1 INTRODUCTION 

The rotational behaviour of pulsars has long been known to be pre- 
dictable, especially in the case of millisecond pulsars (MSPs). Cur- 
rent models suggest that such pulsars have been spun up by ac- 
cretion from their binary companion star to periods of several mil- 
liseconds, making them spin much faster than the more numerous 
younger pulsars, which typically have periods of about a second. 
The rotational stability of MSPs is generally 3-4 orders of magni- 
tude better than that of normal pulsars and on timescales of sev- 
eral years, it has been shown that some MSPs have a timing stabil- 
ity co mparable to the most precise atomic clocks ( Mats akis et al.l 
Il997h . This timing stability is most clearly quantified through 
the technique of pulsar timing, which compares arrival times of 
pulses to a model describing the pulsar, its binary orbit and the 
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interstell ar medium (ISM) be tween the pulsar and Earth (as de- 
tailed bv lEdwardsetal.ll2006l) . This technique has enabled deter- 
mination of physical parameters at outstanding levels of preci- 
sion, such as the orbita l characteristics of binary star systems (e.g. 
Ivan Straten et alj200lh. the masse s of pulsars and their companions 
(e.g. Jacoby et alj|2005l; iNicd 2006) and the turbulent character of 



the ISM (e.g. 



You et al 



1 120071) . The strong gravitational fields of 
pulsars in binary systems have also enabled stringent tests of gen- 
eral relativity (GR) and jdternative theories of gravity, as described 
by, e.g.. lKramer et al.U2006l) and lshat et al.l ( 120081) . Finally, pulsars 
have provided the first evidenc e that gravitational waves (GWs) ex- 
ist at levels predicted by GR jTavlor & Weisber3ll982h and have 
placed the strongest limit yet on the existence of a background of 
GW s in the Galaxy iJenet et al. 2006). 

ISazhirj jl97a) was the first to investigate the potential effect 
of GWs on the times-of-arrival (TOAs) of pulsar pulses and to 
conclude that direct detection of GWs could be possible through 
pulsar timing. Subsequent analyses and theoretical predictions for 
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astronomical sources of GWs have determined that a stochastic 
gravitational- wave background (GWB) from binary black holes in 
the cores of galax ies is the m ost likely signal to be detectable. 
As summarised in Ijenet et 2A\ |_2006), the energy density of such 
a GWB per unit logarithmic frequency interval can be expressed 
as: 



(f) - l^A^l 

■5^0 /rcf 



(1) 



where — 100/ikm s~^ Mpc~^ is the Hubble constant, / is the 
GW frequency, /ref = (1 yr)~^, A is the dimensionless amplitude 
of the GWB and a is the spectral index of the GWB. The one- 
sided power spectrum of the effect of such a GWB on pulsar timing 
residuals is given by: 



P{f) 



-Mff 



where he is the characteristic strain spectrum, defined as: 



hc{f) 



(2) 



(3) 



iJenet et al. also summarized the characteristics and ex- 

pected ranges for various GWBs of interest. Most importantly, for a 
GWB created by supermassive black-hole mergers, a = —2/3 and 
A is predicted to be betwe e n 10~^^ and 10"^"^ jjaffe & Backed 
I2OO3I : IWvithe & Loeb"2003h. ISesana et alj j2008h expanded upon 
these analyses and showed that the actual GWB spectrum strongly 
depends on the merger history, with a variety of spectral indices 
possible. They concluded, however, that a = —2/3 was a rea- 
sonable approximation for practical purposes. For a background of 
GWs that were formed in th e early Universe, a ~ — 1 and the 
amplitude range predicted bv lGrishc hulj ^2005 ) is ^1 = 10~ ^^ — 
10~^^, but standard models (e.g. iBovle & Bu onannol l2008h pre- 
dict much lower amplitudes. A third GWB that may be detected 
by PTAs, is formed by cosmic strings ( Damour & Vilenkin 20q3; 
ICaldwelletalJl996h , with predicted am plitudes between 10 — 
IQ- ^* and a = -7/6 jMaggiordbOOOt) . 

iHellings & DownsI ( 1 19831) first investigated the correlations 



that arise between timing data of different pulsars due to the 
presence of a stochastic and isotropic GWB in the Galaxy. They 
demonstrated that the GWs cause a quadrup olar correlation be- 
tween the timing residuals of different pulsars. [Romanil fl989h and 
[poster & Backe expanded this analysis and introduced the 

concept of a pulsar timing array (PTA), in which an ensemble of 
pulsars is timed and their residuals correlated with each other. The 
PTA co ncept uses the quadrupola r correlation signature first de- 
rived by 'Hellings & DownsI ( Il983h to separate the effect of a GW 
from all other contributions to the residuals, such as intrinsic pulsar 
timing irregularities, clock errors, ISM effects and Solar System 
ephemeris errors. Alternatively, the corr elation signature for non- 
Einsteinian GWs (as recently derived bv lLee et alJbOOSh could be 
used. 

The PTA concept was more rigorously explored bv ljenet et al.l 

( l2005h who first determined the sensitivi ty of PTA experim ents to 
backgrounds of GWs (Equation (12) of Ijenet et al J 120031 Their 
analysis showed that the sensitivity of a PTA depends on four main 
parameters: the number of pulsars, the data span (T), the root- 
mean-square of the timing residuals (simply 'rms' or a henceforth) 
and the number of observations in each of the pulsar timing data 
sets (A'^toa). They further determined that, for a PTA consisting of 
weekly observations of 20 MSPs, all with a timing rms of 100 ns, 
a five-year observational campaign would be required to make a 



~ 3(7 detection of a GWB with a = -2/3 and A = lO"^'^. 
It follows from Equation (12) of Ijenet et all ( |2005[) that the low- 
est amplitude of a GWB from supermassive black-hole mergers to 
which a PTA is sensitive, scales as: 



A^ 



n,GWB OC 



(4) 



Depending on the achievable rms residual of MSPs, an alternative 
PTA could therefore achieve the same results through timing of 20 
MSPs on a biweekly basis for ten years with an rms of close to 
300 ns. This raises two questions related to the potential of PTAs to 
detect a GWB. First, down to which precision can MSPs be timed 
(o"miii) and second, can a low residual rms be maintained over long 
campaigns (i.e. does a/T^''^ decrease with time)? In the context 
of this second question, we will henceforth use the term "timing 
stability" when referring to the potential of an MSP timing data set 
to maintain a constant, preferably low rms residual at all timescales 
up to the timespan of a PTA project, which is typically envisaged 
to be five years or longer. 

It has been shown for a few pulsars that timing with a residual 
rms of a few hundred nanos econds is possible for campaigns last- 
ing a few years. Specificallv. lHotan et alj ( l2006h presented a timing 
rms of 200 ns over two years of timing on PSRs J 17 13-1-0747 and 
J1939-I-2134 (PSR B 1937 +21) and 300 ns ove r two years of tim- 
ing on PSR J1909-3744; ISplaver et all ( l2005h reported a n rms of 
180 ns on six years of timing PSR J17 13-1-0747 and , Verbiest et al.l 
1 I2OO8I) timed PSR J0437-4715 at 200 ns over ten years. Similar 
results for PSRs J0 437-4715 and J1939-I-2134 were obtained by 
iHobbs et al.l ( I2OO8I) over five years of timing. It has, however, not 
been demonstrated thus far that MSPs can be timed with an rms 
residual of 100 ns over five years or more. 

The second question - whether a low rms residual can be 
maintained over ten years or longer, also remains unanswered. 
[Kaspi et al. ( 1 994) detected ex cess low-frequency noise in PSR 
J1939-I-2134; ISplaver et al. 1( 12005 ) presented apparent instabilities 
in long-term timing of PSR J1713-I-0747 and lVerbiest et al.l ( |2008|) 
noted correlations in the timing residuals of PSR J0437— 4715, but 
apart from these, no long-term timing of MSPs has been presented 
to date. Given the low rms residual reported on all three sources, 
it is unclear how strongly the reported non-Gaussian noise would 
affect the use of these pulsars in a GWB detection effort. 

In this article we present the first high-precision long-term 
timing results for a sample of 20 MSPs. The source selection, ob- 
serving systems and data analysis methods are described in fj2] Our 
updated timing models and residual plots for all pulsars in our sam- 
ple are also presented in ij2] allowing the reader a fundamental in- 
spection of the reliability of our timing results. In ^ we perform a 
stability analysis of the timing data, with the dual purpose of iden- 
tifying low-frequency noise in any of our timing data and of as- 
sessing the potential impact of such noise on the use of pulsars in 
a timing array. In Sj4]we outline a new way of quantifying differ- 
ent components of the pulsar timing rms. Through this analysis, 
we separate the levels of receiver noise, noise with a dependency 
on observing frequency and any temporal instabilities, providing a 
bound on the residual rms that might be achievable on a five-year 
timescale. We apply this analysis to three of our most precisely 
timed pulsars. In fj5] we calculate sensitivity curves for ongoing 
and planned PTAs. These sensitvity curves take into account the 
inhomogeneous character of a realistic array (i.e. the fact that the 
rms will differ between pulsars) and assume a bound on residual 
rms as determined in ^ In Sj6]we summarise our findings. 
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2 OBSERVATIONS AND DATA REDUCTION 

2.1 Sample Selection 

The data presented in this paper have been collated from two pul- 
sar timing programmes at the Parkes radio telescope. The old- 
est of the se commenced during the Parkes 70 cm MSP survey 
teailes et al. 1994), aiming to characterise properly the astrometric 
and binary parameters of the MSPs found in the s urvey. Initial tim - 
ing results from thi s campaign were published bv lBell et alJ h997h 
and iToscano et al.l (l999). The bright millisecond pulsars PSRs 
J1713+0747 and B1937+21 (both discovered eariier at Arecibo) 
were also included in this programme. A few years later, as new 
disc overies were made in the Swinburne intermediate-latitude sur- 
vey tedwards et al.ll200]]) . these pulsars were also added, resulting 
in a total of 16 MSPs that were regularly timed by 200 6. Improved 
timin g solutions for these 16 pulsars were presented bv lHotan et al.l 
( [2006 ) and Ord et al. (200©- 

Besides the projects described above, the Parkes Pulsar Tim- 
ing Array (PPTA uManchestej2008l) project commenced more reg- 
ular timing observations of these pulsars in late 2004, expanding the 
number of MSPs to 20 (listed in Table[TJ and adding regular moni- 
toring at a low observing frequency (685 MHz) to allow correction 
for variations of the ISM electron density. A detailed analysis of 
these low frequency observations and ISM effects was recently pre- 
sented bv lYou et alJ t2007) and an a nalysis of the combine d data on 
PSR J0437-4715 was pubhshed bv lVerbiest et al.l j2008l) . For this 
pulsar we will use the timing results presented in that publication; 
for all other pulsars we will present our improved timing models in 

2.2 Observing Systems 

Unless otherwise stated, the data presented were obtained at the 
Parkes 64 m radio telescope, at a wavelength of 20 cm. Two re- 
ceivers were used: the H-OH receiver and the 20 cm multibeam 
receiver jStavelev-Smith et~aL Ifl996). Over the last five years, ob- 
servations at 685 MHz were taken with the 10/50 cm coaxial re- 
ceiver for all pulsars; however, they were only used directly in the 
final timing analysis of PSR 10613—0200, whose profile displays a 
sharp spike at this frequency if coherent dedispersion is applied. For 
PSRs J1045-4509, J1909-3744 and 11939-^2134, the 685 MHz 
observations were used to model and remove the effects of tem- 
poral variations in interstellar dispersion delays and hence are in- 
cluded indirectly in the timing analysis. For all other pulsars any 
such variations were below the level of our sensitivity. 

Three different observing observing systems systems were 
used. Firstly, the Caltech Fast P ulsar Timing Machine (FPTM; 
ISandhuetalJI 19971 ; ISandhJl200lh . between 1994 and November 
2001. This is an autocorrelation spectrometer with a total band- 
width of up to 256 MHz. Secondly, the 256 MHz bandwidth ana- 
logue filterbank (FB) was used in 2002 and 2003. Finally, the 
Caltech-Parkes-Swinburne Recorder 2 (CPSR2; lHotan et al.ll2006l) 
was used from November 2002 onwards. CPSR2 is a baseband data 
recorder with two 64 MHz bandwidth observing bands (one usually 
centred at 1341 MHz, the other at 1405 MHz) and phase-coherent 
dispersion removal occuring in near real time. 

2.3 Arrival Time Determination 

The processing applied differs for data from different observing 
systems. The FPTM data were calibrated using a real-time sys- 



tem to produce either two or four Stokes parameters which were 
later combined into Stokes I. The FB data were produced from a 
search system with no polarimetric calibration possible. This sys- 
tem produced Stokes I profiles after folding 1-bit data. Data from 
both of these systems were integrated in frequency and time to pro- 
duce a single profile for each observation. These observations were 
~25 minutes in duration. For CPSR2 data, in order to minimise 
the effects of aliasing and spectral leakage, 12.5% of each edge of 
the bandpass was removed. To remove the worst radio frequency 
interference, any frequency chaimel with power more than Aa in 
excess of the local median was also removed ("local" was defined 
as the nearest 21 channels and the standard deviation a was deter- 
mined iteratively). CPSR2 also operated a total power monitor on 
microsecond timescales, which removed most impulsive interfer- 
ence. 

The CPSR2 data were next integrated for five minutes and cal- 
ibrated for differential gain and phase to correct for possible asym- 
metries in the receiver hardware. If calibrator observations were 
available (especially in the years directly following the CPSR2 
commissioning, observations of a pulsating noise source, needed 
for polarimetric calibration, were not part of the standard observ- 
ing schedule). Subsequently the data were integrated for the dura- 
tion of the observation, which was typically 32 minutes for PSRs 
12124-3358, J1939-H2134 and I1857-I-0943 and 64 minutes for all 
other pulsars. In the case of PSR 11643— 1224, the integration time 
was 32 minutes until December 2005 and 64 minutes from 2006 
onwards. Finally, the CPSR2 data were integrated in frequency and 
the Stokes parameters were combined into total power. CPSR2 data 
that did not have calibrator observations available were processed 
identically, except for the calibration step. While for some pul- 
sars (like PSR 10437—4715) these uncalibrated data are provably 
of inferior quality (see, e.g. ivan S traten 200^, in our case this is 
largely outweighed by the improved statistics of the larger number 
of TOAs and by the extended timing baseline these observations 
provided. We therefore include both calibrated and uncalibrated ob- 
servations in our data sets. 

Pulse TOAs were determined through cross-correlation of the 
total intensity profiles thus obtained with pulsar and frequency- 
dependent template profiles. These template profiles were created 
through addition of a large number of observations and were phase- 
aligned for both CPSR2 observing bands. As there were only few 
high signal-to-noise observations obtained with the FPTM and FB 
backends for most pulsars, these data were timed against standards 
created with the CPSR2 backend. This may affect the reliability 
of their derived TOA errors. For this reason we have evaluated the 
underestimation of TOA errors for each backend separately, as ex- 
plained in the next section. While the TOA errors were generally 
determined through the standard Fourier phase gradient method, 
the Gaussian interpolation method produced more accurate esti- 
mates for pulsars with low signal-to-noise ratios ( Hotan et al . 200^ 
- specifically for PSRs 10613-0200, 12129-5721, J1732-5049, 
J2124-3358 and J 1045-4509. The PSRCHIVE software package 
dHotan et alj|2004l) was used to perform all of the processing de- 
scribed above. 
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Table 1. Pulsars in our sample. Column 2 gives the reference for the discovery paper, while column 3 provides references to recent or important publications 
on timing of the sources. For the three pulsars with original B1950 names, these names are given beside the J2000.0 names. 



Pulsar 
name 



Discovery 



Previous 
timing 
solution'^ 



Pulse 
period (ms) 



Orbital Dispersion 
period (d) measure 
(cm""^ pc) 



J0437^715 
J06 13-0200 
J07 11-6830 
J1022+1001 
J 1024-07 19 



Johnston et al. (1993) 
Lorimer et al. (1995') 
Bailes et al. (1997) 
Camilo et al. (1996) 
Bailes et al. (1997) 



1,2 
3 

3,4 
3 
3 



5.8 
3.1 
5.5 
16.5 
5.2 



5.7 
1.2 



7.8 



2.6 
38.8 
18.4 
10.3 

6.5 



11045^509 
J1600-3053 
11603-7202 
J 1643- 1224 
J17 13+0747 



Bailes et al. (1994) 
Ord et al. (2006) 
Lorimer et al. (1996) 
Lorimer et al. (1995) 
Foster etal. (1993) 



3 
5 
3 
4 
3,6 



7.5 
3.6 
14.8 
4.6 
4.6 



4.1 
14.3 
6.3 
147.0 
67.8 



58.2 
52.3 
38.0 
62.4 
16.0 



J 1730-2304 
J1732-5049 
J1744-1134 

B1821-24; J1824-2452 
B1855+09; J1857+0943 

J 1909-3744 

B 1937+21; J 1939+2 134 

J2124-3358 

J2 129-5721 

J2145-0750 



Lorimer et al. (1995) 
Edwards & Bailes (2001) 
Bailes et al. (199T) 
Lvneetal. (198T) 
Segelstein et al. (1986) 



Jacobv et al. (2003) 
Backer et al. (1982) 
Bailes etal. (1997) 
Lorimer et al. (1996) 
Bailes etal. (1994) 



4 
7 
3 

8, 10 
3,9 

3, 11 
3,9 
3 
3 

3, 12 



8.1 
5.3 
4.1 
3.1 
5.4 

2.9 
1.6 
4.9 
3.7 
16.1 



5.3 

12.3 
1.5 



6.6 
6.8 



9.6 
56.8 
3.1 
120.5 
13.3 

10.4 
71.0 

4.6 
31.9 

9.0 



° Ref erences: (iVVerbiest etalJ)2008h : (2)lvan Straten et al.'l200lh: (3)lHotan et al.l j2006t) : (4)'Toscano et al.'h 999l): (5)l0rd et al.l j2006l); (6)rSplaver et alj 
b005l): (7) lidw ards & B ailed I2OOII) ; (8) iHobbs et al. (2004); (9) iKaspi et alj h994l) ; (10) Cognard & BackeT l l2004h ; (11) Ijacobv et alj i2005h ; (12) 



iLohmer et al.l )2004l) 
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Table 2. Summary of the timing results, sorted in order of decreasing rms residual. The columns present the pulsar name, the rms timing residual (without 
prewhitening), the length of the data set and the number of TOAs. For PSRs J1939+2134 and J1857+0943 this table only contains the Parkes data. See il2.4l 
and [|3]for details. 



Pulsar 


rms 


T 


Ntoa 


name 


(Ms) 


(yr) 




11 009—3744 


0.166 


5.2 


893 


T171 3+0747 


0.198 


14.0 


392 


J0437-4715 


0.199 


9.9 


2847 


J1744-1134 


0.617 


13.2 


342 


11030+2134 


0.679 


12.5 


168 


Tif^no— '^ns^ 

J 1 UIJU J\U J 


1 12 


6.8 


477 


TOf^l ^— 0900 


1 52 


8 2 


190 


J1824-2452 


1.63 


2.8 


89 


J1022+1001 


1.63 


5.1 


260 


J2145-0750 


1.88 


13.8 


377 


J1643-1224 


1.94 


14.0 


241 


J1603-7202 


1.98 


12.4 


212 


J2129-5721 


2.20 


12.5 


179 


J1730-2304 


2.52 


14.0 


180 


J 1857+0943 


2.92 


3.9 


106 


J1732-5049 


3.23 


6.8 


129 


J07 11-6830 


3.24 


14.2 


227 


J2124-3358 


4.01 


13.8 


416 


J1024-0719 


4.17 


12.1 


269 


J1045-4509 


6.70 


14.1 


401 



PSR J0437-4715 

^1 1 \ ; T" 



PSR J06 13-0200 



PSR J071 1 -6830 






PSR J1022+1001 

O I^T 1 1 p j 1- 



1996 2000 2004 
PSR d1 024-07 19^ 



2000 2004 
PSR t 045-4509 




2008 2000 2008 



2004 




Pm J 1603-7202 

2 FT ■ 



~r 



2007 
>2 

T 1 1 r 1 r- 



±4 



1 — I 1 X — 1, — I — t 1 1_ 



PSf! J 1.643-1 224 




PSR J1713+0747 



PSR J1730-2304 PSR J1 732-5049 

in 






2000 

PSR J 1744-1 134 



20.Q§ 2000 

PSR 41^24-2452 




- un 



T 



2008 2000 tQOa 

PSR J 1857+0943 

■ I I' J I r I 'r J ■! I c t f 1 



2004 2Qm 
PSR J1 909-3744 





- PS 




1996 2000 20D4 fSOS 
PSR J 1939 + 2 134 



2006 2007 Z008 1SS4 1992 2000 mm 
PSR J2 124-335.8 PSR J2 129.-5721 



2004 2007 

PSR J21 45-0750 




1984 1^92. MMM 2008 



2000 



2008 



1996 2Q;Q0 2QQ4 2QQS 



200:0 



2008 



Figure 1. Timing residuals of the 20 pulsars in our sample. Scaling on the x-axis is in years and on the y-axis in /iS. For PSRs J1857+0943 and J1939+2134, 
these plots include the Arecibo data made publically available by Kaspi et al. (1994); all other data are from the Parkes telescope, as described in §2. Sudden 
changes in white noise levels are due to changes in pulsar backend set-up - see §2 for more details. 
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Table 3. Timing parameters for the single pulsars PSRs J071 1-6830, J1024-0719, J1730-2304, J1744-1134, J1824-2452, J1939+2134 and J2124-3358. 
Numbers in brackets give twice the formal standard deviation in the last digit quoted. Note that these parameters are determined using TEMP02, which uses 
the International Celestial Reference System and Barycentric Coordinate Time. As a result this timing model must be modified before being used with an 
observing system that inputs Tempo forniat parameters. See Hobbs et al. (2006) for more information. 



Fit and data set parameters 



Pulsar name 


J0711-6830 


J1024-0719 


J1730-2304 


J1744-1134 


MJD range 


49373.6-54546.4 


50117.5-54544.6 


49421.9-54544.8 


49729.1-54546.9 


Number of TOAs 


227 


269 


180 


342 


rms timing residual (fis) 


3.24 


3.80 


2.52 


0.617 


Reference epoch for P, a. 










5 and DM determination 


49800 


53000 


53300 


53742 


Measured Quantities 


Right ascension, a (J2000.0) 


07:11:54.22579(15) 


10:24:38.68846(3) 


17:30:21.6611(3) 


17:44:29.403209(4) 


Declination, 5 (J2000.0) 


-68:30:47.5989(7) 


-07:19:19.1700(10) 


-23:04:31.29(8) 


-11:34:54.6606(2) 


Proper motion 'ma,fj,a cos 5(mas yr^ ^ ) 


-15.55(8) 


-35.3(2) 


20.27(6) 


18.804(15) 


Proper motion in 5, /is (mas yr^ ^ ) 


14.23(7) 


-48.2(3) 




-9.40(6) 


Annual parallax, tt (mas) 








2.4(2) 


Dispersion measure, DM (cm^'^ pc) 


18.408(4) 


6.486(3) 


9.617(2) 


3.1380(6) 


Pulse frequency, u (Hz) 


182.117234869347(4) 


193.71568356844(13) 


123.110287192301(2) 


245.4261197483027(5) 


Pulse frequency derivative, i/ (10~^® s~-^) . . . 


-4.94406(15) 


-6.95(3) 


-3.05906(10) 


-5.38188(4) 


Prewhitening Terms 


Fundamental wave frequency, ojpw (yr ^ ) . . . 




0.10368 






Amplitude of wave 1 cosine and sine terms. 










^cos,i;Asin,i (10-4 s) 




2(13); 4.7(21) 








Fit and data set parameters 


Pulsar name 


J1824-2452 


J1939+2134 


J2124-3358 




MJD range 


53518.8-54544.9 


46024.8-54526.9 


49489.9-54528.9 




Number of TOAs 


89 


180 


416 




rms timing residual (fis) 


0.986 


0.354 


4.03 




Reference epoch for P, a. 










S and DM determination 


54219 


52601 


53174 




Measured Quantities 


Right ascension, a (J2000.0) 


18:24:32.00796(2) 


19:39:38.561297(2) 


21:24:43.85347(3) 




Declination. 5 (J2000.0) 


-24:52:10.824(6) 


+21:34:59.12950(4) 


-33:58:44.6667(7) 




Proper motion in a, /icx cos S (mas yr^^) .... 




0.072(2) 


-14.12(13) 




Proper motion in 5, /ig (mas yr^ ^ ) 


-9(5) 


-0.415(3) 


-50.34(25) 




Annual parallax, tt (mas) 




0.13(13) 


3.1(11) 




Dispersion measure, DM (cm""^ pc) 


120.502(2) 


71.0227(9) 


4.601(3) 




Pulse frequency, u (Hz) 


327.4055946921(6) 


641.928233642(12) 


202.793893879496(2) 




Pulse frequency deiivative, u {10~^^ s"^) . . . 


-1736.5(3) 


-429.1(6) 


-8.4597(2) 




Prewhitening Terms 


Fundamental wave frequency, ojpw (yr ^ ) . . . 


0.44734 


0.14996 






Amplitudes of cosine and sine terms (IQ— * s): 










wave 1: Acoa.l ; Asin,l 


-20(6); 2.1(14) 


286(41); -413(60) 






wave 2: Acoa,2; ^sin,2 




30(5); 84(12) 






wave 3: Acos.s; ^sin,3 




-21(3); -5.8(9) 






wave 4: Acos,4; Agin, 4 




3.7(5); -2.9(5) 






wave 5: Acos.s; -4sin,5 




0.04(3); 0.68(9) 







Table 4. Timing parameters for binary PSRs J1022+1001, J1600-3053, J1713+0747, J1857+0943, J1909-3744 and J2145-0750. See caption of Table 3 for 
more information. 



Fit and data set parameters 


Pulsar name 


J1022-I-1001 


J1600-3053 


J17 13+0747 


J1857+0943 


J1909-3744 


J2145-0750 


MJD range 


52649.7-54528.5 


52055.7-54544.6 


49421.9-54546.8 


46436.7-54526.9 


52618.4-54528.8 


49517.8-54547.1 


Number of TOAs 


260 


477 


392 


376 


893 


377 


rms timing residual (fis) 


L63 


1.12 


0.198 


1.14 


0.166 


1.88 


Reference epoch for P, a, 5 














and DM determination 


53589 


53283 


54312 


50481 


53631 


53040 


Measured Quantities 


Right ascension, a (J2000.0) 


10:22:58.003(3) 


16:00:51.903798(11) 


17:13:49.532628(2) 


18:57:36.392909(13) 


19:09:47.4366120(8) 


21:45:50.46412(3) 


Declination, S (J2000.0) 


+10:01:52.76(13) 


-30:53:49.3407(5) 


+07:47:37.50165(6) 


+09:43:17.2754(3) 


-37:44:14.38013(3) 


-07:50:18.4399(14) 


Proper motion in a, fi^ cos 6 (mas yr~"^) . 


-17.02(14) 


-1.06(9) 


4.924(10) 


-2.64(3) 


-9.510(7) 


-9.66(15) 


Proper motion in S, /i^ (mas yr~ ^ ) 




-7.1(3) 


-3.85(2) 


-5.46(4) 


-35.859(19) 


-8.9(4) 


Annual parallax, tt (mas) 


1.8(6) 


0.2(3) 


0.94(10) 


1.1(4) 


0.79(4) 


1.6(5) 


Dispersion measure, DM (cm^^ pc) 


10.261(2) 


52.3262(10) 


15.9915(2) 


13.300(4) 


10.3934(2) 


8.9977(14) 


Pulse frequency, v (Hz) 


60.7794479762157(4) 


277.9377070984926(17) 


218.8118404414362(3) 


186.494078620232(2) 


339.31568740949071(10) 


62.2958878569665(6) 


Pulse frequency derivative, u (10~^® s~^) 


-1.6012(2) 


-7.3390(5) 


-4.08379(3) 


-6.20495(6) 


-16.14819(5) 


-1.15588(3) 


Orbital period, Pj, (days) 


7.8051302826(4) 


14.3484577709(13) 


67.825130963(17) 


12.32719(4) 


1.533449474590(6) 


6.83893(2) 


Orbital period derivative, Pb (10~^^) 






41(20) 


3(3) 


5.5(3) 


4(3) 


Epoch of periastron passage, Tq (MJD) . . . 


53587.3140(6) 


53281.191(4) 


54303.6328(7) 


50476.095(8) 




53042.431(3) 


Projected semi-major axis, x = asin j (s) 


16.7654074(4) 


8.801652(10) 


32.3424236(3) 


9.230780(5) 


1.89799106(7) 


10.1641080(3) 


x(10~i4) 


1.5(10) 


-0.4(4) 






-0.05(4) 


-0.3(3) 


Longitude of periastron, luq (deg) 


97.75(3) 


181.85(10) 


176.190(4) 


276.5(2) 




200.63(18) 


Orbital eccentricity, e(10~^) 


9.700(4) 


17.369(4) 


7.4940(3) 


2.170(6) 




1.930(6) 


K = esina;o (10~*) 










-0.4(4) 




rj = e cos u)o (10~®) 










-13(2) 




Ascending node passage, Tasc (MJD) .... 










53630.723214894(4) 




Periastron advance, ui (deg/yr) 












0.06(6) 


Sine of inclination angle, sin i 


O.73I 


0.8(4) 




0.9990(7) 


0.9980(2) 




Inclination angle, i (deg) 


472- 




78.6(17) 








Companion mass, Mc (Mq) 


1.05^ 


0.6(15) 


0.20(2) 


0.27(3) 


0.212(4) 




Longitude of ascending node, f2 (deg) .... 






67(17) 









Table 5. Timing parameters for binary PSRs J0613-0200, J1045-4509, J1603-7202, J1643-1224, J1732-5049 and J2129-5721. See caption of Table 3 
for more information. 



Fit and data set parameters 


Pulsar name 


J0613-0200 


J1045-4509 


J1603-7202 


J1643-1224 


J1732-5049 


J2129-5721 


MJD range 


51526.6-54527.3 


49405.5-54544.5 


50026.1-54544.7 


49421.8-54544.7 


52056.8-54544.8 


49987.4-54547.1 


Number of TOAs 


190 


401 


212 


241 


129 


179 


rms timing residual (/xs) 


1.52 


6.70 


1.98 


1.94 


3.23 


2.20 


Reference epoch for P, a, S 














and DM determination 


53114 


53050 


53024 


49524 


53300 


54000 


Measured Quantities 


Right ascension, a (J2000.0) 


06:13:43.975142(11) 


10:45:50.18951(5) 


16:03:35.67980(4) 


16:43:38.15544(8) 


17:32:47.76686(4) 


21:29:22.76533(5) 


Dechnation, S (J2000.0) 


-02:00:47.1737(4) 


-45:09:54.1427(5) 


-72:02:32.6985(3) 


-12:24:58.735(5) 


-50:49:00.1576(11) 


-57:21:14.1981(4) 


Proper motion in a, /ic< cos <5 (mas yr~ ) . 


1.84(8) 


-6.0(2) 


-2.52(6) 


5.99(10) 




9.35(10) 


Proper motion in 5, i^s (tnas yr~^) 


-10.6(2) 


5.3(2) 


-7.42(9) 


4.1(4) 


-9.3(7) 


-9.47(10) 


Armual parallax, tt (mas) 


0.8(7) 


3.3(38) 




2.2(7) 




1.9(17) 


Dispersion measure, DM (cm""' pc) 


38.782(4) 


58.137(6) 


38.060(2) 


62.409(2) 


56.822(6) 


31.853(4) 


Pulse frequency, u (Hz) 


326.600562190182(4) 


133.793149594456(2) 


67.3765811408911(5) 


216.373337551614(7) 


188.233512265437(3) 


268.359227423608(3) 


Pulse frequency derivative, u (10"^*^ s""^) 


-10.2308(8) 


-3.1613(3) 


-0.70952(5) 


-8.6438(2) 


-5.0338(12) 


-15.0179(2) 


Orbital period, (days) 


1.1985125753(1) 


4.0835292547(9) 


6.3086296703(7) 


147.01739776(5) 


5.262997206(13) 


6.625493093(1) 


Epoch of periastron passage. To (MJD) . . . 


53113.98(2) 


53048.98(2) 




49577.9689(13) 




53997.52(3) 


Projected semi-major axis, x = a sin i (s) 


1.0914444(3) 


3.0151325(10) 


6.8806610(4) 


25.072614(2) 


3.9828705(9) 


3.5005674(7) 


i(10-") 






1.8(5) 


-4.9(5) 




1.1(6) 


Longitude of periastron, uo (deg) 


54(6) 


242.7(16) 




321.850(3) 




196.3(15) 


Orbital eccentricity, e (10~^) 


0.55(6) 


2.37(7) 




50.578(4) 




1.21(3) 


K — esincjo (10~®) 






1.61(14) 




2.20(5) 




77 — ecosuJo (10~^) 






-9.41(13) 




-8.4(4) 




Ascending node passage, Tgac (MJD) 






53309.3307830(1) 




51396.366124(2) 
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2.4 Timing Results 

The TEMP02 software package jHobbs etal .l2006h was used to cal- 
culate the residuals from the TOAs and initial timing solutions (Ta- 
blelTJ. In order to account for the unknown instrumental delays and 
pulsar-dependent differences in observing setup, arbitrary phase- 
offsets were introduced between the data from different backends. 
Where available, data at an observing frequency of 685 MHz were 
included in an initial fit to inspect visually for the presence of dis- 
persion measure (DM) variations. In the case of PSRs J1045— 4509, 
J1909-3744, J1939-I-2134 and J0437-4715, such variations were 
significant and dealt with in t he timing s oftwar e through a method 
similar to that presented by I You et The average DM 

values presented in Tables |3] |4] and |5] were determined from the 
20 cm data exclusively. The uncertainties in these DM values do not 
take into account possible pulse shape differences between the pro- 
files at these slightly varying frequencies. We updated all the pul- 
sar ephemerides to use International Atomic Time (implemented 
as TT(TAI) in TEMP02) and the DE405 Solar System ephemeris 
( IStandisril2004h . 

In order to correct for any underestimation of TOA uncer- 
tainties resulting from (amongst others) the application of CPSR2- 
based template profiles to the FB and FPTM data (as mentioned 
in W1.2i\ and to allow comparison of our timing model parameters 
to those published elsewhere, the TOA uncertainties were multi- 
plied by error factors (so-called "EFACs") that are dependent on 
the pulsar and observing system. Specifically, this part of the anal- 
ysis was performed as follows. First the timing data from each ob- 
serving system were prewhitened by fitting harmonically related 
sine/cosine pairs if required. Next the TOA uncertainties were mul- 
tiplied by an EFAC value that produced a reduced value of unity 
for that prewhitened subset of the data. Because of potential non- 
Gaussian noise in the data, application of these backend-specific 
EFACs does not necessarily result in a reduced value of unity for 
the entire, recombined data set. To account for such non-Gaussian 
noise in the data, a 'global' EFAC was applied to the entire data set, 
making the reduced after prewhitening equal to unity and in- 
creasing the parameter uncertainties reported in the timing models 
appropriately. As mentioned the prewhitening method was based 
on fitting of sine/cosi ne pairs to the data, according to t he following 
formula described bv lMartiiJ l200lh (and replicated in lHobbs et al.l 
l2006h : 

A.R = ^ Asin.k sin(fctjp„ At) + A^os.k cos(fctjp„Af) 
fc=i 

where k runs over all sine/cosine pairs, nn is the total number of 
harmonically related pairs fitted, ^cos.fc and A^i^.k are respectively 
the amplitude of the fe* cosine and sine waves and Wpw is the fun- 
damental frequency derived from: 

_ 27r 
'^•'"^ r(l + 4/nH) 

with T the length of the data set. If prewhitening terms were in- 
cluded in the final fit, we provide the values for ojpw, Acos,k and 
^sin.fc as part of our timing model. 

Because the potential non-Gaussian noise present in the data 
is the subject of our investigations in the remaining sections of this 
paper, the prewhitening terms as well as the global EFACs were not 
included in any subsequent analysis. The residuals plotted in Figure 
[T]and the parameters presented in Table[2]therefore do not include 
prewhitening terms or global EFACs. 

The system-specific EFACs were generally less than two, with 



the only major outliers being the CPSR2 data of PSR J 1939-1-2 134 
with an EFAC of 5.27 and the 32-min CPSR2 integrations (pre- 
2004 CPSR2 data) of PSR J 1643 -1224, which have an EFAC 
of 4.9. In the former case this large EFAC may be caused by in- 
complete prewhitening, as the non-Gaussian noise is badly mod- 
elled by polynomials or sine/cosine pairs. The underestimation 
of PSR J 1643— 1224 TOA uncertainties is likely caused by the 
low signal-to-noise ratio of these observations, which causes the 
Fourier phase gradient method to under estimate TOA errors (as 
previously reported bv lHotan et alj|200^ . We note that the EFAC 
for the 64 minute integrations is much lower, at 2.5. In deriving the 
timing models, the global EFAC was at most 1.1 and for more than 
half of our sources less than 1.05. 

The fact that most of our EFAC values are close to unity and 
show little variation with backend, suggests that the parameter and 
error estimates are fairly robust. In order to account for the different 
sensitivity of the backends used and to limit effects of scintillation 
on our timing, we opt for a weighted analysis. It is therefore impor- 
tant to consider the impact of the TOA errors and applied EFACs 
on the different parts of this analysis. Given that for most pulsars 
the EFACs applied to the different backends are nearly equal, the 
resulting timing models will be little affected by these EFAC val- 
ues. The reported uncertainties on the timing model parameters will 
be affected but will be comparable to previous publications, since 
ou r analysis method is similar. A full error analysis (as suggested 
by I Verbiest et al.P 2008') is needed to provide any more reliable pa- 
rameter uncertainties. Since the focus of this present paper is on 
the overall timing stability and implications for pulsar timing ar- 
ray science, we defer such error analysis (and the interpretation of 
any previously unpublished parameters in our timing models) to a 
later paper. We have, however, investigated the effect of weight- 
ing and EFACs on the timing stability analysis (^, but have not 
uncovered any unexpected deviations beyond statistical noise. We 
therefore conclude that the weighting and applied EFACs do not 
invalidate our analysis. 

A summary of the lengths of the data sets and the achieved rms 
residual can be found in Table (2] highlighting the superior resid- 
ual rms of PSRs J1909-3744, J0437-4715 and J1713-I-0747 when 
compared to other pulsars. The timing residuals for our data sets are 
presented in Figure[T]and the timing models are presented in Tables 
[3]|4]and|5] where 2a errors are given, in accordance with previous 
practice. We encourage observers to use the improved models when 
observing. We also note that all but a few of the parameters in our 
timing models are consistent with those published previously. 



3 PULSAR TIMING STABILITY 

In SjT] we demonstrated that one of two vital questions relating to 
the potential of PTAs to detect a GWB is whether a low residual 
rms can be maintained over long timespans (a property we refer to 
as "timing stability"). Effectively, this question breaks down into 
two parts: to what degree of significance low-frequency noise is 
present in our pulsar timing data and how any such low-frequency 
noise can be expected to affect sensitivity to a GWB. In order to an- 
swer this question fully, a spectral-analysis-based investigation of 
pulsar timing residuals that includes identification and modelling 
of potential non-Gaussian noise sources, would be required. Be- 
cause of various pulsar timing-specific issues such as clustering of 
data, large gaps in data sampling and large variations in error-bar 
size, however, standard spectral analysis methods fail to provide 
reliable power spectra of pulsar timing data. We therefore use the 
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alternative app r oach p rovided by the statistic, as described by 
iMatsakis et alj ( Il997h . A brief explanation of this statistic, along 
with a presentation of the CTz values of our data is presented in i)3.1l 
and a discussion of these results in terms of PTA-science is pro- 
vided in 



steeper than a simulated GWB, implying that this pulsar will most 
likely not be very useful for long-term PTA projects, although its 
low rms residual on short time spans might make it useful for de- 
tection of burst-type sources. 



3.1 (Tz Stability Analysis 

Originally proposed bv lMatsakis et alj ( ll997h . the CTz statistic is de- 



fined as: 



^2 



' 2\l/2 



2^5 

where () denotes the average over subsets of the data, C3 is deter- 
mined from a fit of the polynomial 

Co +C-i{t- to) + C2{t - to)'^ + Cz(t - tof 

to the timing residuals for each subset and r is the length of the 
subsets of the data. In order for the CTz values to be independent 
of each other, we use r = T, r/2, T/A, T/8, . . . only. The inter- 
pretation of this statistic i n terms of power spe ctra deserves some 
attention. As presented bv lMatsakis et al.l ( ll997l) . a power spectrum 
with spectral index 

P{u) cc f 
would translate into a curve: 

(Tz(r) oc r'', 

where the spectral indices are related as: 

,= (-^(^ + 3) (5) 
1—2 otherwise. 

Equation ^ implies that spectra have different slopes in a CTz 
graph than in a power spectrum. Along with the CTz graphs for our 
data sets, Figure[2]provides some examples of spectra for guidance: 
lines with a slope of —3/2 (dotted lines in Figure|2l represent spec- 
trally white data (/3 = into Equation l|5j gives /i — —3/2) and 
a GWB with a spectral index a = —2/3 in the gravitational strain 
spectrum (and therefore a spectral slope /3 = —13/3 in the timing 
residual spectrum, as follows from equations ^ and l|3j) would 
have a positive slope of 2/3 in CTz (dashed lines). 

Comparison of such theoretical slopes to the data is, however, 
non-trivial since the data are strongly affected by effects from sam- 
pling and fitting. As an illustration of such effects, the top-left plot 
of Figure |2] shows two ctz curves derived from simulations. The 
first one is the full line that approximates the ctz curve for PSR 
J1713+0747. In this case the ctz values of 1000 simulations of white 
noise with the timing rms and sampling of the PSR J1713+0747 
data set were averaged. This curve is not perfectly parallel to the 
theoretical curve with slope —3/2 due to sampling, varying TOA 
uncertainties and model fitting. Comparison of the white noise sim- 
ulations with the actual PSR JI7I3+0747 data indicates that there is 
not a significant, steep red-noise process affecting the timing resid- 
uals for this pulsar. The second simulation in the top-left plot of 
Figure |2] is the dot-dashed line, which is the average ctz graph of 
2000 simulations of white noise with an artificial GWB and the 
sampling of the PSR JI939+2134 data set, fitted for pulse period 
and spindown. This simulated curve does not reach the theoretical 
slope of 2/3 because of flattening off at low frequencies caused by 
sampling, fitting and leakage resulting from these. This simulation 
also demonstrates that the PSR J 1939+2 1 34 curve is significantly 
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Figure 2. stability parameter for tlie 20 pulsars in our sample. The dotted slanted lines represent white noise levels of 100ns (bottom) and lO^s (top); 
the dashed slanted line shows the steepness introduced to pre-fit residuals by a hypothetical GWB (see JT); pulsars whose curve is steeper than this line (like 
PSR J1939+2134), can therefore be expected to be of little use to PTA efforts on long timescales. The top left figure further shows the average ctz values 
resulting from 1000 simulations of white noise residuals sampled at the times of the PSR J17 13+0747 data set and fitted for the PSR J17 13+0747 timing 
model parameters (full line). This demonstrates that the PSR J1713+0747 data do not - within the sensitivity provided by the (Tz statistic - contain a significant, 
steep red-noise process. The dash-dotted line in the top left figure shows the average of 2000 simulations for white noise combined with a GWB, sampled at 
the times of the PSR J 1939+2 134 data set and fitted for pulse period and period derivative. These simulated results provide an example of the combined effect 
sampling and model fitting can have on the (Jz statistic, even in the case of white noise. 
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Figure 3. CTz graphs for all pulsars in our sample, showing both weighted (black and green lines) and unweighted (red and blue lines) results, as well as results 
including EFACs (black and red lines) and excluding EFACs (green and blue lines). The dashed lines represent theoretical white noise at levels of 100 ns and 
10 fis. The EFAC values of the PSR J0437— 4715 data were lost in processing, so the green and blue curves are missing for that particular pulsar For all other 
pulsars all four curves are present, though they do frequently overlap. 
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The (Tz graphs of our data are shown in Figure (2] A compar- 
ison of these curves to those obtained from an unweighted analy- 
sis or from an analysis that does not contain the EFAC values de- 
scribed in i]2.4l is presented in Figure|3] a colour version of which is 
available online. This graph demonstrates that use of weighting or 
EFACs does not affect the data in any statistically significant man- 
ner, other than to decrease the effect of the white noise component 
in case of scintillating pulsars. 

Comparison of the PSR J1713+0747 and PSR J1939+2134 
data with the simulated curves shown, along with the invariability 
of the shape of CTz plots to weighting or application of error factors, 
shows that the CTz parameter provides a good first-order discrimina- 
tion between pulsars that do not exhibit significant, steep red noise 
(like PSR J1713+0747) and those that do have timing instabilities 
which could mask a GWB (like PSR J 1939+2 134). 

3.2 Timing Stability Conclusions 

Figure|2]shows that PSR J 1939+2 134 has red noise with a level and 
steepness that will limit its use in GWB-detection efforts that last 
more than about two years. Four other pulsars (PSRs J0613— 0200, 
J1024-0719, J1045-4509 and J1824-2452) show some indica- 
tion of similar red noise, but longer timing and lower white-noise 
levels are needed to determine this with statistical significance. For 
all the other pulsars we have no evidence that the red noise that may 
be present in the timing residuals below the white-noise level has 
a spectral index that prevents GWB detection on timescales of five 
years to a decade. 

We have been unable to detect timing instabilities with an am- 
plitude and spectral slope that could mask a GWB in the timing 
data of PSRs J1713+0747 and J1744-1134, notwithstanding their 
long data spans and low timing rms which should make them highly 
sensitive to any low-frequency noise. Using Equation l|4j, it can be 
shown that the data sets of PSRs J1713+0747 and J0437-4715 al- 
ready meet the requirements for a ten-year long PTA experiment, 
proving that at least for some pulsars the timing stability and rms 
residual required to detect a GWB on timescales of ten years or 
more is achievable. The challenge for such long-term projects will 
therefore be to find more pulsars like these, or to replicate these 
results for other existing pulsars, by increasing the sensitivity of 
observing systems. 

It must be noted that the study of irregularities in pulsar tim- 
ing data (often referred to as "timing noise") can be much more 
extensive than presented here. Given our main aim of assessing the 
impact on PTA science and the absence of clear timing noise in 
most of our data sets, precise modelling or bounding of timing ir- 
regularities as well as a thorough discussion of the potential sources 
of any observed timing instabilities, has not been included in this 
analysis. 



4 ANALYSIS OF RESIDUAL RMS 

As an alternative to the long-term PTA detection efforts discussed 
in the previous section, a shorter-term detection is possible if tim- 
ing at lower residual rms is achievable. The standard scenario for a 
(relatively) short-term GWB detection by a PTA requires five years 
of weekly observ ations with a timing rms of 100 ns for 20 MSPs 
jjenet et al.ll20Cl3) . Since a residual rms of 100 ns has never been 
maintained over five years, the possibility that some intrinsic prop- 
erty of MSPs induces instabilities at that level, remains open. In 
this section, we will address that issue by evaluating how much the 



timing rms of some of our most precisely timed pulsars may be 
reduced. 

We separate three different categories of contributions to the 
pulsar timing residuals: 

Radiometer noise o-Rad: The Gaussian noise component that 
scales with the radiometer equation and which is mainly deter- 
mined by the shape and signal-to-noise (S/N) ratio of the observed 
pulse profiles. 

Frequency-systematic effects o^' This category of noise contri- 
butions contains most effects that produce timing residuals depen- 
dent on the observing frequency. This includes interstellar effects 
such as interstellar scintillation and DM variations. 

Temporal-systematic effects o-ri This category contains all 
time-dependent effects such as calibration errors, instabilities in 
the observing systems, clock errors, errors in the Solar System 
ephemerides, GWs and intrinsic pulsar timing noise. 

As it is impossible to get direct measures of the three contri- 
butions listed above, we base our analysis on the following three 
measurements: 

Total timing rms o: This is simply the timing residual rms of 
the data considered. It contains all three effects: 

O"^ = o-Rad + + O-?- (6) 

Sub-band rms a^h- In ^A.2\ we will introduce this new measure 
which is l/-\/2 times the weighted rms of the offset between the 
residuals of two simultaneous observing bands with different centre 
frequencies. Since the observations in the two observing bands are 
simultaneous, their offset is determined by the radiometer noise and 
by frequency-systematic effects (as the observing bands are centred 
at slightly different frequencies). We can therefore write: 

o^sb = cr|ad + o-'i. (7) 

Theoretical radiometer noise crpad : In i]4.1| we will calculate 
(jRad directly from the pulse profiles used in our timing. 

Using these three measures and equations (|6} and Q, the three 
contributions to the timing residuals can be isolated, the results of 
which are described in i]4.3l 

Our analysis will be based on the CPSR2 data of PSRs 
J1909-3744, J1713+0747 and 11939+2134. We restrict this analy- 
sis to the CPSR2 data because it is of superior quality to the data of 
older backend systems (see ^2.2\ and because it consists of the five 
most densely sampled years of observations. We focus on three of 
the most precisely timed pulsars in order to obtain the best limits 
on achievable residual rms. In doing so, we omit PSR J0437— 4715 
because the advanced calibration schemes u sed in its analysis (see 
IVerbiest et ai]|2008l : Ivari"stratenll2004l I2OO6) complicates our ef- 
forts and beca use reported non-Ga ussian noise in the timing data 
of this pulsar ("Verbie st et al.ll200^ may imply an inferior limit to 
that derived from PSRs 11909-3744 and J1713+0747. Note that 
the purpose of this analysis is to uncover the potential limit for 
high-precision timing: it is already known (see e.g. ^ that MSPs 
have different amounts of time-dependent noise, so the limit we 
will derive from PSRs 11909-3744 and J1713+0747 does not have 
to hold for all MSPs. However, it does suggest that other pulsars 
may achieve similar rms residual and that a PTA-size sample of 20 
MSPs at such rms residual may mainly depend on increased sensi- 
tivity of current observing systems and new discoveries in ongoing 
and future surveys. 
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Table 6. Breakdown of weighted timing residuals for three selected pulsars. Given are the total timing rms of the ~5 years of CPSR2 data (cr), the sub-band 
timing rms (o-gb), the radiometer noise (a-Had)i the temporal systematic ((Tt) and the frequency systematic (cr^) contributions to the timing rms. All values are 
in ns and apply to 64min integrations. See [|4]for more information. 



Pulsar name 


cr 


0"sb 


""Had 


0"t 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


J1909-3744 


166 


144 


131 


83 


60 


J17 13+0747 


170 


149 


105 


82 


106 


J1939+2134 


283 


124 


64 


254 


106 



4.1 Theoretical Estimation of Radiometer Noise, o-Rad 

The level at which the radiometer noise adds to the timing residu- 
als can be determined based on the pul sar's observed pul se profile 
shape and brightness, as described by Ivan Stratei] j20065 . Equa- 
tion (13) of that publication provides the following measure (notice 
we only consider the total intensity, So, to allow direct comparison 
with our timing results): 



CTRad = P X = P : 



E 



(8) 



where is the m frequency of the Fourier transform of the 
pulse profile, is the total power at that frequency, i;o is the 
white noise variance of the profile under consideration, A'^ is the 
total number of time bins across the profile and A'max is the fre- 
quency bin where the Fourier transform of the pulse profile reaches 
the white noise level, <;o. V i& the expected variance in the phase- 
offset or residual, P is the pulse period and o-Rad is the residual rms 
predicted for the pulse profile considered. 

In order to use Equation {Sj on our data, we first integrated all 
our pulse profiles together, weighted by signal-to-noise ratio, after 
which Equation ^ was applied to the final profile. Subsequently 
ORad was renormalised to 64min integrations through use of the 
radiometer equation. In order to check this result, we also applied 
the equation to all individual pulse profiles contained in this anal- 
ysis and averaged the results in a weighted way - resulting in the 
same answer, which is given in column four of Table [6] The value 
for PSR J1909— 3744 shows that even at this low residual rms, ra- 
diometer noise dominates the timing rms. Applying this method to 
the other MSPs in our sample, we found that almost all our timing 
residuals are dominated by radiometer noise. For more than half 
of our sample of 20 MSPs, URad is of the order of a microsecond 
or more. This demonstrates the need for longer integration times, 
larger bandwidth and/or larger collecting area. 

4.2 Estimating Frequency-Dependent Effects 



As described in ii2.2l the CPSR2 pulsar backend records two ad- 
jacent, 64 MHz-wide frequency bands simultaneously. This allows 
determination of a unique measure of a sub-set of timing irregular- 
ities, which we will call the "sub-band rms", a^h'- 



71^ E.i/e: 



2 

z.mn 



(9) 



where the sums run over all observing epochs i, r^.m and ri,n are 
the residuals of either observing band (named m and n respectively) 

at epoch i and ei^mn = yc^^^^ + e^ is the average TOA error at 

epoch i. Effectively, the sub-band rms is l/\/2 times the weighted 



rms of the offset between the residuals of the two bands. This im- 
plies it contains all contributions to the total rms that are not time- 
dependent but either statistically white or dependent on the observ- 
ing frequency, as described earlier. Note, however, that many of 
these effects have both a temporal and frequency component. Given 
our sampling, it should therefore be understood that (specifically in 
the case of DM variations) only part of these effects is contained in 
(J I,, while the remaining contributions are contained in Ot- 

The sub-band rms for the three selected pulsars is presented in 
column three of Table[6] 

4.3 Discussion 

Based on equations ^ and and the three measures cr, a^h and 
ORad determined in the preceding paragraphs, the three contribu- 
tions to the rms (uRad, <^v and CTt) can now be estimated. Their 
values are presented in columns 4, 5 and 6 of Table |6] In order to 
assess the potential for 100 ns timing of these sources over a five- 
year timescale, we will now discuss the possible means of reducing 
these three contributions. 

The radiometer noise o-Rad scales for different telescopes or 
observing systems according to the radiometer equation: 



o-Rad OC 



A^s\/Bt 



(10) 



where B is the bandwidth of the observing system used, t is the 
integration time, A^s = »7^^-j- is the effective collecting area of 
the telescope (with -q the aperture efficiency and D the telescope 
diameter) and Tsys is the system temperature of the receiver. 

The frequency systematic contributions are not as easily 
scaled for different observing systems, but t hey can be de- 
creased and research on thi s front is progressing jYou et alj|2007l: 
iHemberger & Stinebrin j2008 :.Walker et ah. 20081). Also, by reduc- 
ing the radiometer noise, any measurements of DM variations will 
become more precise, which will enhance corrections for these ef- 
fects and therefore decrease the contribution of Ou. We also note 
that since these effects are frequency dependent, the employment of 
very large bandwidth receivers or coaxial receiver systems such as 
the 10/50 cm receiver at the Parkes observatory, may lead to highly 
precise determination and correction of these effects. Furthermore, 
increased collecting area and bandwidth may enable future tim- 
ing observations at higher observational frequencies, which would 
limit the size of these effects. We therefore suggest that ai, does 
not ultimately limit the achievable rms residual, but may largely 
be corrected for if current research and technological development 
progress. 

The wide variety of sources that add to the temporal system- 
atic make predictions about its future evolution hard. Sources such 
as intrinsic pulsar timing noise are (as yet) impossible to mitigate. 
Errors in the terrestrial clocks or in the Solar System ephemerides 
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are expected to decrease as better models become available or 
as timing arrays provide their own improved solutions for these 
models. Instabilities in the observing system m ay to some degree 
be m itigated by improved calibration methods jvan Stratenll2004 
l2006h . Simultaneous observations of a single source at multiple 
observatories may also lead to detection and correction of instru- 
mental instabilities and the time-dependent effect of DM variations 
may also be mitigated, as explained above. 

Following from the above, we stress the fact that all contribu- 
tions to and (jRad may be mitigated, but that certain contribu- 
tions to (Tr cannot be corrected. This implies that this last class of 
effects will ultimately limit the residual rms that can be reached. 
We will therefore use the temporal systematic contribution to the 
rms ((Jt, column five in Table |6j as an upper limit on the poten- 
tial rms residual of the MSPs under investigation. Note that this 
is a conservative upper limit since significant protions of Or may 
be expected to be mitigated. However, without relative quantifica- 
tion of the various contributions to Ot, this limit cannot be reliably 
decreased. 

Given the discussion above, we note that the potential timing 
residual rms of PSRs J1909-3744 and J1713+0747 is predicted to 
be below 100 ns on a five-year timescale. This implies that the stan- 
dard scenario of 100 ns timing over five years is possible provided 
techniques currently being developed for mitigation of frequency- 
dependent effects are successful, more sensitive observing systems 
are used and more bright, stable MSPs like PSRs J1909— 3744 and 
J17 13+0747 are found. 



5 PROSPECTS FOR GRAVITATIONAL WAVE 
DETECTION 

Ijenet et alj ^20051) derived the expected sensitivity of a PTA to 
a GWB with given amplitude, A, both for homogeneous arrays 
(where all pulsars have comparable timing residuals) and inhomo- 
geneous arrays. They also pointed out the importance of prewhiten- 
in£|the residuals to increase sensitivity at larger GWB amplitudes. 
In the current section, we will build upon their analysis to provide 
more realistic predictions for ongoing and future timing arrays. We 
extend their analysis in three fundamental ways. 

Firstly, we use the rms timing residuals presented in Table [2] 
These results provide an inhomogeneous set of rms's with a realis- 
tic spread. We assume the residuals are statistically white and will 
therefore not change with the timescale of the timing array project. 
Our analysis in Sj3] shows that for most pulsars this assumption is 
reasonable, especially on timescales of order five years. 

Secondly, w e do not apply exactly the same algorithm as 
Ijenet et al] j2005l) . In Appendix [A] we present a derivation of PTA 
sensitivity to a GWB in a manner that provides some guidance on 
analysing the data. We assume that the prewhitening and correla- 
tion are handled together by computing cross-spectra and we esti- 
mate the amplitude of the GWB directly rather than using the nor- 
malised cross correlation function. We assume that the non-GWB 
noise is white, but ca n be different for e ach pulsar. Our results are 
very close to those of Ijenet et and using our method we 

successfully reproduced the scaling law. Equation The analysis 



^ In this context, prewhitening refers to a technique that flattens the power 
spectrum of a time series by means of weighting. This flattening optimises 
the sensitivity of a PTA to steep red spectra such as those introduced by a 
GWB. 



could be easily extended to include non-white noise if a model for 
the noise were available. 

Finally, in order to generalise the results from our Parkes data 
to telescopes in other parts of the world, we scale the residuals 
based on realistic parameters for various PTA efforts listed in Ta- 
ble [T] In doing so, we scale cTRad (see according to Equation 
dlOl ). As discussed in i]4.3l some improvements in and ar can 
be expected in coming years, especially as the radiometer noise is 
decreased. While quantification of any such improvement is practi- 
cally impossible, we will apply the same radiometer scaling to cr^ 
as we apply to o-Rad and assume Or to be constant at 80 ns for all 
pulsars at all telescopes. This may provide a slight disadvantage 
for larger telescopes, but overall we consider this a reasonable yet 
conservative approach. 

5.1 Ongoing PTA Projects 

We consider five ongoing PTA efforts: 

Current: Refers to the data presented in this paper, using the 
longest overlapping time span of the sample: five years. This ig- 
nores the fact that the PSR 11824-2452 data set is shorter, but this 
globular cluster pulsar may not prove useful in a PTA project last- 
ing longer than a few years anyway. We therefore assume that a 
replacement is found and has identical timing rms over a time span 
of five years. 

Predicted PPTA: Assumes the usage of 256 MHz of bandwidth 
at the Parkes telescope, which implies a four-fold bandwidth in- 
crease and therefore a two-fold decrease in timing rms. The PPTA 
is the only one to be considered for more than five years, mainly 
in order to demonstrate the large impact a doubling of campaign 
length can have, but also because several years of high precision 
timing data with that bandwidth do already exist ( lManchestej2008l) 
for all 20 MSPs. 

NANOGrav: Assumes Arecibo gain for the ten least well-timed 
pulsars and GBT gain for the ten best-timed pulsars, in order to get 
a fairly equal rms for all 20 MSPs. (Since we consider an upper 
limit on the rms residual, the advantage of Arecibo over the GBT 
is limited for the brightest sources.) 

EPTA: Assumes monthly observations with five lOOm-class 
telescopes dlanssen et al.ll2008h . 

EPTA-LEAI0: Interferometrically combines the five tele- 
scopes of the EPTA to form a single, larger one. This decreases 
the number of observations, but increases the gain. 

An important caveat to this analysis is that several of the pul- 
sars under consideration cannot be observed with most Northern 
telescopes, because of the telescope declination limits. We there- 
fore assume stable MSPs to be discovered in the Northern hemi- 
sphere. As mentioned before, we also assume that progress will be 
made in the mitigation of frequency-dependent ISM and calibration 
effects. Finally, this analysis is based on the Parkes data presented 
in this paper and therefore assumes systematic effects to be at most 
at the level of the Parkes observing system used. 

The sensitivity curves presented in Figure |4] seem to justify 
cautious optimism for GWB detection through PTA experiments 
on timescales of five to ten years, provided current models of 
GWBs are correct. While none of the curves in Figure |4] reach 
the minimum predicted GWB Amplitude of 10^^^ at a detection- 
significance level of three, their sensitivity can be expected to in- 

^ Large European Array for Pulsars 
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crease up to an order of magnitude through extension of the cam- 
paigns to a decade-long timescale, as illustrated by the difference 
between the "Predicted PPTA" and "Current sensitivity" curves. 
The GWB predictions may, however, change if other effects such 
as black-hole binary stalling occur. The models do, furthermore, 
rely on a substantial number of poorly determined input parame- 
ters, such as what fraction of galaxy growth happens by merging 
dSesana et al.] |2009'). Since only the merging of galaxies results in 
binary black holes and hence contributes to the GWB, this mass 
fraction is crucial for any reliable prediction of GWB strength. 

As explained in i]4.3l the temporal systematic contribution to 
the rms, ct^, is a conservative upper limit to the ultimate residual 
rms. In this analysis of PTA efforts, however, we have used the 
value of 80 ns as a hard lower limit on the timing rms, a. This limits 
the potential for reduction of the rms and explains the equivalence 
of the NANOGrav and EPTA-LEAP efforts. Finally, the strong de- 
pendence on the timescale, T, of the project underscores the impor- 
tance of timing stability analysis over much longer time spans and 
continued observing. While our ctz analysis on PSR J1713-I-0747 
provides the first evidence for high timing stability over timescales 
beyond ten years, such timing stability must still be demonstrated 
for many more MSPs. 



5.2 Future PTA Projects 

With the completion of the Square Kilometre Array (SKA) 
pathfinders expected in three years time, we consider the poten- 
tial of the Australian SKA Pathfinder (ASKAP), the South African 
Karoo Array Telescope (MeerKAT) and the Chinese Five hundred 
meter Aperture Spherical Telescope (FAST) for PTA programmes. 
ASKAP is primarily designed for Hi surveys and therefore sac- 
rifices point source sensitivity for a wide field of view, whereas 
MeerKAT's design is better suited for point source sensitivity over 
a more limited field of view. FAST is an Arecibo-type single dish 
with a total diameter of 500 m of which 300 m is illuminated, re- 
sulting in a substantially larger sky coverage than is possible with 
Arecibo. The expected architecture for these telescopes is listed in 
Table|7]- notice we assume phase-coherent combination of the sig- 
nals of all ASKAP and MeerKAT dishes, effectively resulting in 
a single telescope of diameter 107 m for MeerKAT and 76 m for 
ASKAP. 

The resulting sensitivity curves are drawn in Figure |5] along 
with a hypothetical curve for the most sensitive telescope currently 
operational, the Arecibo radio telescope. This figure clearly shows 
the advantage MeerKAT holds over ASKAP for PTA work, in num- 
ber of dishes, bandwidth and system temperature. The sensitivity 
of Arecibo is much higher than that of either interferometric pro- 
totype and is just inferior to FAST. As for the NANOGrav and 
EPTA-LEAP projects analysed earlier, the advantage of FAST over 
Arecibo is strongly limited by the bound of 80 ns we imposed on 
the achievable rms residual. 

Note the usefulness of Arecibo for PTA work is limited by 
the restricted sky coverage and hence available pulsars. While both 
MeerKAT and ASKAP can see large parts of the Southern sky, the 
sky coverage of Arecibo as well as the short transit time make an 
exclusively Arecibo-based PTA practically impossible; however, its 
potential as part of a combined effort (Figure|4l( or in a global PTA, 
is undeniable if the level of systematic errors is small compared to 
the radiometer noise. As for any Northern telescope, the usefulness 
of FAST will mostly depend on the discovery of good timing MSPs 
at positive declinations. 



6 CONCLUSIONS 

We have presented the first long-term timing results for the 20 
MSPs constituting the Parkes Pulsar Timing Array (PPTA). We 
have shown that only PSR J1939-I-2134 has timing instabilities that 
limit its use for long-term GWB efforts, while the PSR J17 13+0747 
data already demonstrate the requirements for GWB detection on 
a timescale of ten to fifteen years are achievable. Overall, the tim- 
ing stability of the investigated MSPs was found to be encouraging 
even though potential timing instabilities were detected in four pul- 
sars (in addition to PSR J1939+2134). 

It was demonstrated that even on our most precisely timed 
MSPs, white noise is a dominant contribution, suggesting that our 
residual rms will be much improved with current wide-bandwidth 
systems. We placed a conservative upper limit of -^80 ns on intrin- 
sic timing instabilities that will ultimately limit the residual rms. 
We interpreted this result in the context of ongoing and future pul- 
sar timing array projects, demonstrating the realistic potential for 
GWB detection through pulsar timing within five years to a decade, 
provided technical and data reduction developments evolve as ex- 
pected. For PTA efforts in the Northern hemisphere, the discov- 
ery of bright and stable MSPs in the Northern sky will be cru- 
cial. Given the location of currently known MSPs, the prospects 
of the MeerKAT SKA pathfinder as a gravitational-wave detector 
are found to be particularly good. 



APPENDIX A: PTA SENSITIVITY 

In this Appendix we derive a simplified formalism for estimating 
the sensitivity of a pulsar timing array (PTA) to a stochastic and 
isotropic gravitational-wave background (GWB) of given ampli- 
tude, A. This derivation produce s results eq uivalent to those result- 
ing from equation ( 14) of Jenet et al.l j2005l) . but is more readily im- 
plemented and inherently treats optimal weighting (or prewhiten- 
ing) of the pulsar power spectra. 

The detection statistic is the sample cross-covariance of the 
residuals of two pulsars i and j, separated by an angle 9i,j: 

1 ^ 

(where ri(t) is the residual of pulsar i at time t, Ns is the num- 
ber of samples in the cross covariance and T is the data span). 
The expected value of R{9ij) is the covariance of the clock error, 
which is 100% correlated, plus the cross covariance of the GWB, 
''"gwC(^ij ) - The clock error can be included in the fit, but one must 
also include its variance in the variance of the detection statistic. It 
is better to estimate the clock error and remove it, which also re- 
moves its "self noise". So in subsequent analysis we neglect clock 
noise and effects of possible Solar System ephemeris errors. We 
model the pulsar timing residuals as a GWB term and a noise term: 
r{t) — rGw(t) + rtq{t), with varia nces a% and a%. C( Si.i) is the 
cross-correlation curve predicted bv lHellings & DownsI ( Il983h . as 
a function of the angle between the pulsars, 6ij: 

C{^^,d) ^ ^^^°g^^ I + ^ 

in which a:: = (1 — cos 0i.j)/2. 

Since the detection significance will be limited by the variance 
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Figure 4. Sensitivity curves for different PTA efforts. Note tlie "NANOGrav" and "EPTA - LEAP" curves are practically coincident. Gravitational waves are 
predicted to exist in the range 10^15 - 10-1"'. See text and Table|2]for more information. 

Table 7. Assumed parameters for future and ongoing PTA efforts. Besides the names of the different PTAs, the columns contain the number of telescopes 
A^Tcl. the observing bandwidth B, the telescope diameter D, aperture efficiency r\, system temperature Tgys, observing regularity and the duration of the 
project, T. 



PTA 


A^Tol 


B 


D 






Observing 


T 


name 




(MHz) 


(m) 




(K) 


regularity 


(yrs) 


Current 


1 


64 


64 


0.6 


25 


weekly 


5 


Predicted PPTA 


1 


256 


64 


0.6 


25 


weekly 


10 


NANOGrav 


2 


256 


305; 100 


0.5; 0.7 


30; 20 


monthly 


5 


EPTA 


5 


128 


100 


0.7 


30 


monthly 


5 


EPTA - LEAP 


11 


128 


224 


0.7 


30 


monthly 


5 


Arecibo 


1 


512 


305 


0.5 


30 


two-weekly 


5 


FAST_^ 


1 


400 


500 


0.36 


20 


two-weekly 


5 


ASKAP^ 


40 


256 


12 


0.8 


50 


weekly 


5 


MeerKAT_^ 


80 


512 


12 


0.7 


30 


weekly 


5 



" Under the LEAP initiative, five lOOm-class telescopes will be combined into an effective 224 m single telescope. 
[Nan 1 2006); Jin et al. 1 2008) 
http://www.atnf.csiro.au/projects/askap/specs.html 
|jonas|j2007l) 

in the sample cross covariance, we consider After prewhitening this becomes (notice our notation crpw = Q)'- 
Var(i?(6l,,j)) 

= Var ((rcw.z + rN.«)(rGw,, + rN.j)/iVs)) Var(7?pw (61,,^)) 

2 2 (1 + C(^»,j)^) , '^N.iO-G,j + 0-G,i'^N,j , CrN,iO-N,j 4 ( 1 + C(^i.j )^) , 2 (t'N.i + ^N.^) , .On 

— '^G,iO"G,j Ti ^ ^7 1 Ti • ^ 6:G ^ £>g ^ 



iVs A^s iVs Ns Ns 



(A3) 



(A2) 
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Figure 5. Sensitivity curves for the two main SKA pathfinders, Arecibo and FAST. Gravitational waves are predicted to exist in the range 10"^^ - 10"". 
See discussion in f|5]and TabIe[2]for more information. 



in which we have used qq ^ — qq ^ — Qq, which will be proved 
shortly. 

We derive the gravitational-wave power from equations {3} 
and for a GWB with spectral index a = —2/3: 



(A4) 



with K a constant proportional to the amplitude of the GWB and 

/rcf = 1 J/r"\ 

Defining the corner frequency, /c, as the frequency at which 
the gravitational wave power equals the noise power, enables us 
to use equation ( IA4t to determine the noise power: PNoisc = 

As illustrated bv ljenet et alj ( 120051) , the steep spectral index 
of GWB-induced residuals implies that large gains in sensitivity 
can be achieved through optimal prewhitening of the data. Assess- 
ment of the variance of both the GWB and noise components of 
the residuals after prewhitening, can most easily be done through 
integration of the spectral powers, multiplied by the whitening fil- 
ter, W{f), which is a type of Wiener filter, designed to minimize 
the error in the estimation of ao and is of the form: W{f) = 
Pgwb / (Pgwb + PnoIsc)^ - Reseating the weighting function thus 
defined, we get: 



W{f) = C- 



-13/3 



(l + (///c)-^3/3)^ 



(A5) 



with C a normalisation constant chosen for convenience to be: 

-26/3 



(///ref)- 



r (1 + (///c)-"/')^ 



(A6) 



The prewhitened variances then become: 



2 

Qg 



2 

Qn 



-13/3 



/ 

K 



4" 1 



(A7) 



(A8) 



which justifies our choice for C and shows that, based on our 
weighting scheme, Qc^i — Qa^j — K, as used earlier. 

Since the spectra are effectively bandlimited to /c after 
prewhitening, both the GWB and noise will have the same number 
of degrees of freedom, namely: A'dof = iTohsfc — 1, where Tobs 
is the length of the data span and therefore the inverse of the low- 
est frequency, implying there are Tobs/c independent frequencies 
measured below /c. Since each frequency adds a real and imagi- 
nary part, there are twice as many degrees of freedom as there are 
independent frequency samples; quadratic fitting removes a single 
degree of freedom from the total. Notice that \/ Ndof,iNdoi,j is the 
number of independent samples in the cross-covariance spectrum 
and therefore replaces A's in equations iAH and l |A3l l. 

The optimal least-squares estimator for K (and hence for the 
amplitude of the GWB), based on a given set Rpw{di.j) with un- 
equal errors, is (from equations JAU and ((AT}) : 



K : 



Efipw(6>i,j)C(gi,j)/Var(J^Pw,»,j) 
EC(e>j)VVar(i?Pw,»,,) 



(A9) 
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The variance of this estimator is: 



Var(A') 



(AlO) 



EC(»ij)VVar(i?Pw,ij) 

We can now write the expected signal-to-noise of a given tim- 
ing array as the square root of the sum over all pulsar pairs of equa- 
tion \A1\ divided by the square root of equation l lAlOt 



iVpsr-l iVpi 

E E 



£'GC^\/^dof,iA'"dof,j 



Rewriting leads to: 



s = 



E E 



3 

(All) 



(A12) 



where q'^ = gyi,i/QG- 
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